Projecting the development of antimicrobial resistance in Neisseria gonorrhoeae from antimicrobial surveillance data: a mathematical modelling study

Background The World Health Organization recommends changing the first-line antimicrobial treatment for gonorrhoea when ≥ 5% of Neisseria gonorrhoeae cases fail treatment or are resistant. Susceptibility to ceftriaxone, the last remaining treatment option has been decreasing in many countries. We used antimicrobial resistance surveillance data and developed mathematical models to project the time to reach the 5% threshold for resistance to first-line antimicrobials used for N. gonorrhoeae. Methods We used data from the Gonococcal Resistance to Antimicrobials Surveillance Programme (GRASP) in England and Wales from 2000–2018 about minimum inhibitory concentrations (MIC) for ciprofloxacin, azithromycin, cefixime and ceftriaxone and antimicrobial treatment in two groups, heterosexual men and women (HMW) and men who have sex with men (MSM). We developed two susceptible-infected-susceptible models to fit these data and produce projections of the proportion of resistance until 2030. The single-step model represents the situation in which a single mutation results in antimicrobial resistance. In the multi-step model, the sequential accumulation of resistance mutations is reflected by changes in the MIC distribution. Results The single-step model described resistance to ciprofloxacin well. Both single-step and multi-step models could describe azithromycin and cefixime resistance, with projected resistance levels higher with the multi-step than the single step model. For ceftriaxone, with very few observed cases of full resistance, the multi-step model was needed to describe long-term dynamics of resistance. Extrapolating from the observed upward drift in MIC values, the multi-step model projected ≥ 5% resistance to ceftriaxone could be reached by 2030, based on treatment pressure alone. Ceftriaxone resistance was projected to rise to 13.2% (95% credible interval [CrI]: 0.7–44.8%) among HMW and 19.6% (95%CrI: 2.6–54.4%) among MSM by 2030. Conclusions New first-line antimicrobials for gonorrhoea treatment are needed. In the meantime, public health authorities should strengthen surveillance for AMR in N. gonorrhoeae and implement strategies for continued antimicrobial stewardship. Our models show the utility of long-term representative surveillance of gonococcal antimicrobial susceptibility data and can be adapted for use in, and for comparison with, other countries. Supplementary Information The online version contains supplementary material available at 10.1186/s12879-023-08200-4.

Susceptible individuals become infected with either type of infection k ∈ {1, 2} according to a force of infection βSI k . Infected individuals recover spontaneously at rate ν, or through treatment. Here we introduce prescription data under the form of the forcing function p(t), representing the probability that a prescription at time t actually includes the antibiotic of interest. If the antibiotic of interest is prescribed (probability p(t)) then the effect of treatment differs according to the resistance status. Infections by wild-type bacteria (compartment I 1 ) can either recover (with probability 1 − µ) or develop resistance in one single step (with probability µ). For infections by resistant-type bacteria (compartment I 2 ) recovery through treatment occurs with a lower efficacy ϵ ∈ [0, 1]. If another antibiotic is prescribed (probability 1 − p(t)), recovery through treatment occurs at the same rate τ for compartments I 1 and I 2 .

Reparameterization of β
We assume an initial situation where the prevalence of N. gonorrhoeae is at endemic equilibrium. To do that, we reparameterize β as a function of the other parameters, and of the prevalence at endemic equilibrium I * .
In a first step, we retrieve the formula of R 0 in the single-step model using the next generation matrix method (Diekmann et al, 1990;Van den Driessche and Watmough, 2002). This implies the computation of the infection matrix F : and of the migration matrix V : We can then find R 0 as the largest eigenvalue of the next generation matrix F V −1 , which results in: In a second step, we consider that the prevalence at endemic equilibrium I * is related to R 0 through: We can now express the transmission parameter β as function of I * and the other parameters: We thus replace β by this formula in the model, ensuring that the model will start at endemic equilibrium. From now, we treat I * as a parameter, with a normal prior distribution based on estimates of prevalence ranging between 0.16% and 0.38% in HMW and between 1.19% and 2.79% in MSM (Fingerhuth et al. 2016).

Joint model and inference
In total, this model describes the population-level dynamics of resistance development in relation to antibiotic usage in a population with five parameters: {I * , ν, τ, µ, ϵ}. In addition, we consider the initial conditions of resistance in the population, that we treat as another parameter ρ.
To support identifiability, we expand this framework to jointly model the development of resistance in HMW (i = 1) and MSM (i = 2). This takes advantage of the fact that some parameters may be assume to have common values in these two groups, such as the recovery rate ν, the probability of mutation given treatment µ and the treatment efficacy ϵ. The other parameters are left independent (τ i , I * i and rho i ) between the two groups, except for the constraint that τ must be higher for MSM than for HMW (implemented using the ordered data type in Stan).
We estimate all parameters by fitting the model to resistance data using following likelihood: where k t,i is the number of isolates with resistance at time t in population i, n t,i is the sample size, and κ is an overdispersion parameter.

Resistance
We argue that using a prior on µ that favors smaller values allows for a wider range of scenarios, and is thus more adequate that the "flat" priors.

Implementation
We implemented this model in Stan, and conduct parameter inference with Hamiltonian Monte Carlo using Stan default NUTS algorithm. We assessed the quality of the inference by applying diagnosis tools (divergent transitions, tree depth, E-BFMI), and by observing the trace plots and the posterior predictive check plots.
The model code is available in models/binary_model.stan. The R function to format the data and actually run inference is available in fit_binary_model.R. The application to GRASP data is done in main-binary.R.

Model description
We consider an alternative multi-step model that treats resistance acquisition as a multi-step, cumulative process. The model follows a similar structure as the single-step model but includes multiple infected compartments {I 1 , . . . , I k } instead of two. These K compartments represent increasing levels of antibiotic resistance and correspond to the K MIC classes reported in GRASP (e.g., K = 8 in the case of ceftriaxone): This model relies upon two central assumptions. First, the probability of developing one more step of resistance upon treatment µ is the same for every class. Second, increasing levels of AMR lead to a linear decrease in treatment efficacy, with a linear interpolation between ϵ 1 (fixed to 100%) and ϵ K ∈ [0, 1] (estimated): A progressive decrease of treatment efficacy with MIC is compatible with the pharmacodynamical concept of a "period with the free drug level above MIC' ' necessary to achieve treatment efficacy [?]. This second model is also based on five parameters: {β, ν, τ, µ, ϵ K }.

Joint model and inference
We use the same reparameterization of β in terms of the other parameters and the prevalence at endemic equilibrium I * . The initial conditions of resistance are now modelled by a simplex vector of K elements ρ. We also jointly model the growth of MIC in HMW (i = 1) and MSM (i = 2), with the same common parameters ν, µ and ϵ K .
We estimate all parameters by fitting the model to resistance data using following likelihood: where k t,i is the number of isolates within each MIC class at time t and in population i, and ϕ is an overdispersion parameter.

Prior distributions
We selected the following weakly-informative prior distributions: The slightly larger prior on µ is justified by the fact that the probability of developing one more step of resistance upon treatment has to be larger than the probability of developing full resistance directly. We validate our choice of priors with prior predictive checks: load("models/samples_2022-05-17/S_multistep_grasp_azithro_2022-05-17.Rdata") plot_summary_prior2(SIM_multistep_grasp_azithro,lim=2050,colmic = "Greys") Year

Resistance
We see that the chosen priors allow for a large range of scenarios.

Implementation
We also implemented this model in Stan and use the same inference procedure. The model code is available in models/multistep_model.stan. The R function to format the data and actually run inference is available in fit_multistep_model.R. The application to GRASP data is done in main-multistep.R.

Removing data from 2009-2010
We considered a sensitivity analysis where data from 2009 and 2010 was removed. The objective was to evaluate the impact of the temporary rise of the number of observed cases of high resistance in these years, that can be attributed to the international circulation of a multi-drug resistance clone. This was simply implemented by removing observations from 2009 and 2012 from the data.

Increasing prevalence
Another sensitivity analysis aimed at relaxing our hypotheses regarding the prevalence of N. gonorrhoeae infections. In the main analysis, we assumed an initial situation where the prevalence of N. gonorrhoeae is at endemic equilibrium, and reparameterized β as a function of the other parameters. We now consider a situation where the prevalence of N. gonorrhoeae is steadily rising over time. This is implemented by adding a time-dependent component to β, so that it starts at the pre-specified prevalence, and rises linearly throughout the years: where t is the number of years since initiation. This is implemented in models binary_modelB.stan and multistep_modelB.stan. We considered three scenarios with the parameter ζ fixed to 0.001, 0.005 or 0.01, corresponding to increasing slopes in the rise of prevalence.